module fft_wrap

!----------------------------------------------
! VARIABLES
!----------------------------------------------

  use dimensions, only : jx, jy, jz, jv

  integer(kind=8), parameter :: FFTW_ESTIMATE = 0
  integer(kind=8), parameter :: FFTW_MEASURE  = 1
  integer(kind=8), parameter :: FFTW_FORWARD = -1
  integer(kind=8), parameter :: FFTW_BACKWARD = 1

  integer(kind=8) ::  plan_xforward,  plan_xbackward
  integer(kind=8) ::  plan_yforward,  plan_ybackward
  integer(kind=8) ::  plan_zforward,  plan_zbackward
  integer(kind=8) ::  plan_vforward,  plan_vbackward

  complex ::  xIn(jx),  yIn(jy), zIn(jz), vIn(jv)
  complex ::  xOut(jx),  yOut(jy), zOut(jz), vOut(jv)

!----------------------------------------------
! SUBROUTINES
!----------------------------------------------


contains


!-- fftw initialization --

subroutine fft_init()

  use mpi_wrap

  implicit none

  write(msg,*) char(10), ' Initializing fft_wrap module'
  call mpi_wrap_msg

  call dfftw_plan_dft_1d(plan_xforward, jx, xIn, xOut, FFTW_FORWARD, FFTW_ESTIMATE)
  call dfftw_plan_dft_1d(plan_yforward, jy, yIn, yOut, FFTW_FORWARD, FFTW_ESTIMATE)
  call dfftw_plan_dft_1d(plan_zforward, jz, zIn, zOut, FFTW_FORWARD, FFTW_ESTIMATE)
  call dfftw_plan_dft_1d(plan_vforward, jv, vIn, vOut, FFTW_FORWARD, FFTW_ESTIMATE)

  call dfftw_plan_dft_1d(plan_xbackward, jx, xIn, xOut, FFTW_BACKWARD, FFTW_ESTIMATE)
  call dfftw_plan_dft_1d(plan_ybackward, jy, yIn, yOut, FFTW_BACKWARD, FFTW_ESTIMATE)
  call dfftw_plan_dft_1d(plan_zbackward, jz, zIn, zOut, FFTW_BACKWARD, FFTW_ESTIMATE)
  call dfftw_plan_dft_1d(plan_vbackward, jv, vIn, vOut, FFTW_BACKWARD, FFTW_ESTIMATE)

end subroutine


!-- forward transform --

subroutine fft_forward(n, datain, dataout)

  implicit none

  integer  n
  complex  datain(:)
  complex  dataout(:)


  if (n .eq. jx) then

     call dfftw_execute_dft(plan_xforward, datain, xOut)
     dataout = xOut

  else if (n .eq. jy) then

     call dfftw_execute_dft(plan_yforward, datain, yOut)
     dataout = yOut

  else if (n .eq. jz) then

     call dfftw_execute_dft(plan_zforward, datain, zOut)
     dataout = zOut

  else if (n .eq. jv) then

     call dfftw_execute_dft(plan_vforward, datain, vOut)
     dataout = vOut

  else

     write (*,*) "wrong size of data in fft_forward"

  end if

end subroutine



!-- forward transform --

subroutine fft_backward(n, datain, dataout)

  implicit none

  integer  n
  complex  datain(:)
  complex  dataout(:)


  if (n .eq. jx) then

     call dfftw_execute_dft(plan_xbackward, datain, xOut)
     dataout = xOut

  else if (n .eq. jy) then

     call dfftw_execute_dft(plan_ybackward, datain, yOut)
     dataout = yOut

  else if (n .eq. jz) then

     call dfftw_execute_dft(plan_zbackward, datain, zOut)
     dataout = zOut

  else if (n .eq. jv) then

     call dfftw_execute_dft(plan_vbackward, datain, vOut)
     dataout = vOut

  else

     write (*,*) "wrong size of data in fft_backward"

  end if


end subroutine

end module
